NT@UW-i3-o8, INT-PUB-13-004 



On the use of the Discrete Variable Representation Basis in Nuclear Physics 
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The discrete variable representation (dvr) basis is nearly optimal for numerically representing wave 
functions in nuclear physics: Suitable problems enjoy exponential convergence, yet the Hamiltonian 
remains sparse. We show that one can often use smaller basis sets than with the traditional harmonic 
oscillator basis, and still benefit from the simple analytic properties of the dvr bases which requires 
no overlap integrals, simply permit using various Jacobi coordinates, and admit straightforward 
analyses of the ultraviolet (uv) and infrared (ir) convergence properties. 
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F"|roblems in nuclear physics typically require solv- 
* ing the one-body Schrodinger equation in three- 
dimensions. Numerically representing wavef unctions 
requires limiting both ultraviolet (uv) and infrared (ir) 
scales: a finite spatial resolution (i.e., a lattice) character- 
izes the highest representable momenta A, while a finite 
size (i.e. a cubic box of volume L 3 ) determines the largest 
physical extent. Nuclear structure calculations are his- 
torically dominated by the use of the harmonic oscillator 
(ho) basis of ho wave functions. The appeal of the ho 
basis stems from the shape of the self-consistent field 
obtained for small nuclei, which can be approximated by 
a harmonic potential at small distances from the center 
of the nucleus. One can also use the Talmi-Moshinsky 
transformation to separate out the center-of-mass mo- 
tion in products of single particle ho wavefunctions. 
Recent efforts have been made to determine a minimal 
ho basis set, and to understand its convergence and 
accuracy [ij|2j|. 

Here we advocate that the discrete variable representa- 
tion (dvr) - in particular the Fourier plane-wave basis - 
enjoy most of the advantages of the ho basis, but with a 
significant improvement in terms of computational effi- 
ciency and simplicity, thereby admitting straightforward 
uv and ir convergence analyses and implementation. 

Consider wavefunctions in a cubic box of volume 
L 3 with momenta less than A. The total number of 
quantum states in such a representation is given by the 
following intuitive formula - the ratio of the total phase 
space volume to the phase space volume of a single 3D 
quantum state: 
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One obtains the same result [3] using Fourier analysis: 
there are exactly Nqs linearly independent functions in 
a cubic 3D box of volume L 3 with periodic boundary 
conditions and wave-vectors less than k c = A/h in each 
direction. These can be conveniently represented in the 
coordinate representation with N equally spaced points 
in each direction and lattice constant a = 7t/k c = 7th/A = 



L/N for a total of N 3 = Xq S coefficients. The maximum 
wave-vector k c is simply the Nyquist frequency [3I - one 
gains nothing by sampling the functions on intervals 
("times") finer than a. 

The wavefunctions can also be represented in mo- 
mentum space using a discrete fast Fourier transform 
(fft) I4]. The momentum representation consists of Nqs 
coefficients on a 3D cubic lattice with spacing 27th/L and 
extent —A < y x ,y,z < A. Using the fft to calculate spa- 
tial derivatives is not only fast, but extremely accurate - 
often faster and more accurate than finite difference for- 
mulas. We use an even number of lattice points (N = 2 n 
is best for the fft) and quantize the three momenta 

(Px,y,z = hk. X/ y /Z ) 
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The Fourier basis uses plane waves - e.g. exp(iknx) in 
the x-direction - but these can be linearly combined to 
form an equivalent sinc-function basis: 



i|) k (x) =sinck c (x-x k ) 
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This is similar to the difference between Bloch and Wan- 
nier wave functions in condensed matter physics. An ad- 
vantage of this basis is that it is quasi-local 4>k(xi) = 5^ 
allowing one to represent external potentials as a diago- 
nal matrix V^i ~ V(xi c )5] C t (19) . 

The plane wave basis can thus be interpreted as a peri- 
odic dvr basis set, which has been discussed extensively 
in literature (see |3]-(8l and the references therein), and 
one can take advantage of Fourier techniques and the 
useful dvr properties. 

In general, dvr bases are characterized by two scales: 
a uv scale A = hk c defining the largest momentum 
representable in the basis, and an ir scale L defining 
the maximum extent of the system. In many cases, the 
basis is constructed by projecting Dirac delta-functions 
onto the finite-momentum subspace: For example, the 
sinc-function basis Q is precisely the set of projected 
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Dirac delta functions 4>nW = Pp^A 
subspace |p| ^ A [5 -7]. (It can be non-trivial, however, 
to choose a consistent set of abscissa maintaining the 
quasi-locality property.) The basis thus optimally covers 
the region [—L/2, L/2) x [—A, A) for each axis in phase 
space, and leads to an efficient discretization scheme 
with exponential convergence properties. 

The dvr basis admits a straightforward analysis of 
the uv and ir limits, allowing one to construct effective 
extrapolations to the continuum and thermodynamic 
limits respectively. The uv effects may be analyzed 
by simply considering the properties of the projection 
Pp^A used to define the basis, and the ir limit for peri- 
odic basis is well understood by techniques like those 
derived by Beth, Uhlenbeck, and Liischer ||9][ioj|. We 
would like to emphasize an additional technique here: 
The ir limit is characterized by 27tft/L - the smallest 
interval in momentum space resolvable with the ba- 
sis set. For some problems, one can efficiently circum- 
vent this limitation by using "twisted" boundary condi- 
tions n|)(r + L) = exp(i0B)4>(r) or Bloch waves as they are 
known in condensed matter physics. In particular, av- 
eraging over 0g € [0, 27t ) will completely remove any ir 
limitations (without changing the basis size) for periodic 
and homogeneous problems, effectively "filling-in" the 
momentum states p n ^ p n + b-0B/L < Pn+i • Extensions 
of these formulas to the case of a box with unequal sides 
is straightforward. 

To demonstrate the properties of the dvr basis, we 
contrast it with the ho basis. The periodic dvr ba- 
sis (plane-waves) shares the ease of separating out the 
center-of-mass. In particular, one can use Jacobi coor- 
dinates to separate out the center-of-mass motion with- 
out evaluating Talmi-Moshinsky coefficients, leading to 
simpler and more transparent implementations. The 
quasi-locality of the dvr basis offers an additional im- 
plementation advantage over the ho basis: one need 
not compute wavefunction overlaps to form the poten- 
tial energy matrix. In contrast with the ho basis, the 
kinetic energy matrix K is no-longer diagonal, but it has 
an explicit formula 1 24 1, and is quite sparse, unlike the 
potential energy operator in the ho basis. 

Consider the ho wavefunctions with energy L ^ 
ha)(N +3/2): the maximum radius and momenta are 
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where b = A/ft/mtu is the oscillator length. For large N, 
N w RA/2ru Thus, to expand a wavefunction with extent 
2R containing momenta |p| < A requires at least 
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states. To contrast, the dvr basis covering the required 
volume of phase space Q with L = 2R and A is 
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Since these states are localized, one can further im- 
pose Dirichlet boundary conditions, allowing functions 
only of the type sin(k n x) with knL = nn (instead of 
exp(ik n x)), thereby keeping only half of the momenta: 
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Choosing a cubic box with Dirichlet boundary con- 
ditions, sides L = 40 fm, and maximum momentum 
A = 30oMeV/c gives 
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a somewhat surprisingly small number of states. For 
symmetric states, one could further the reduce the basis 
by imposing cubic symmetry, decreasing the basis size 
by another factor of 8. 

Finally, one can fully utilize spherical symmetry with 
a related Bessel-function dvr basis gaining a factor of 
tt/6, and thereby besting the ho basis 



0.8 < 1. 
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In this counting, spin and isospin degrees of freedom 
which occur in both bases have been omitted. 

The Bessel-function dvr basis set |5f|7} El follows 
from a similar procedure of projecting Dirac delta func- 
tions for the radial Schrodinger equation. The angular 
coordinates are treated in the usual manner using spheri- 
cal harmonics, but the radial wavefunctions are based on 
the Bessel functions (see Refs. [7I ITi ] for details) which 
satisfy the orthogonality conditions 
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where z voc = k c r va - the zeros of the Bessel functions 
Iviz^a.) = - define the radial abscissa r V/0l . The dvr 
basis set is 
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Integral and differential operators have simple forms 
in the dvr basis (see Refs. [5-7] and the code Bi2l[i3| ). 
In principle, a different basis (and corresponding ab- 
scissa) should be used for each angular momentum 
quantum number v; In practice, good numerical accu- 
racy is obtained using the v = basis jo(zon r /R) and 
the v = 1 basis )-\ (zi n r/R) respectively for even and odd 
partial waves ETT1IT21 . In the s-waves case, the abscissa 
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are simply the zeros of the spherical Bessel function 

j (z) =sin(z)/z: 



ZOr 



n.7i, 



T17T 



n = 1,2,3,. ..,N, (13) 



and correspond to the id basis with Dirichlet boundary 
conditions mentioned earlier. The zeros for j 1 (z) lie 
between the zeros of jo(z). The number of dvr functions 
needed to represent with exponential accuracy a radial 
wavefunction is 
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to be compared (in the limit N — > 00) with 

Rk c 
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(15) 



In the last formula we have divided by an additional 
factor of 2, since N = 2n + I changes in steps of 2. 

A major drawback of the ho wavefunctions that is 
rarely mentioned is that, for modest values of N and I 7^ 
0, the radial wave functions concentrate in two distinct 
regions: around the inner and outer turning points of 
the effective potential V(r) = ri 2 l(l + 1 )/2mr 2 +ma> 2 r 2 /2. 
By adding components with larger values of N, one 
modifies the wavefunction at both small and large dis- 
tances, leading to slow convergence. In contrast, the dvr 
functions are concentrated around a single lattice site. 
Thus, adding more components only affects the solution 
in the vicinity of the additional lattice points leaving the 
states largely unaffected elsewhere. 

For nuclei one can gain insight with some estimates. 
Cutoffs of A = 6ooMeV/c and R = 1 .5 • • • 2A 1 / 3 fm 
should satisfy most of the practical requirements, lead- 
ing to 



0.7- -0.8A 1/6 fm, 



(16a) 
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60- --80 A - 1/3 MeV, (16b) 



compared to the value 40 A -1 / 3 MeV one finds in typical 
monographs J14I. Using only half the value of A = 
30oMeV/c naturally halves the value of Tia). 

We end with demonstrations of the dvr method [jTJj. 
We start with the harmonic oscillator problem in id 
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4>(x) = E(Kx), (17) 



where we choose the harmonic oscillator frequency 
according to Eq. < |i6b[ , varying the lattice constant 
a = 7r/k c and L = Nq. The dvr method is sometimes 
referred to as the Lagrange method in numerical anal- 
ysis |8[, and functions are usually represented on the 
spatial lattice 
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Potential matrix elements usually have a simple and 
unexpectedly accurate representation (quasi-locality) 
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(19) 



where the functions fkM are a linear combination of 
plane-waves and form an orthonormal set (these formu- 
lae apply for even numbers of abscissa as required by 
efficient implementations of the fft) 
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where x k and p n were defined in Eq. 1(2}. As before, the 
functions fk(xi) are simply the normalized projections of 
the periodic Dirac functions on the dvr subspace [5H7J' 
and satisfy 



Y_ afk(*n)fl(Xn) = §kl- 



(22) 



The sinc-function basis Q is obtained in the limit N — > 
00 (if a = 1). Derivatives can be computed using the 
projected derivative of the Dirac delta function g x (x), 
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and is exact for all momenta except for the highest one 
Pmax = hn/a. Some care must be taken when dealing 
with this state, but physics should not depend on this if 
A is properly chosen. 

While the potential matrix is diagonal, the dvr kinetic 
energy is a matrix in coordinate representation: 



Kiel 



mN 1 ? sin 2 Tzjk-i) 



6ma 2 



1 + 



k = l. 



(24) 



This matrix is full matrix in id, but sparse in 3D where 
only 1 /N 2 of the matrix elements are non-vanishing. The 
ho Hamiltonian 1 (17} is thus represented in the dvr basis 
with periodic boundary conditions as 
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The implementation of Dirichlet boundary conditions 
uses the v = Bessel function basis (see the matlab 
code [12] for 1 = and also Ref. [8 J for other possible 
dvr basis sets in id). 
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Figure 1. (color online) Difference in spectrum between the 
dvr Hamiltonian p^j l and the ho energies (n + 1/2)htu. The 
three (blue) curves with pluses have fixed uv scale (lattice 
constant a = 1, k c = n/a) with L =e {30,40,50} and cu = 2n/L 
from left to right. The (red) curves with dots have fixed L = 30 
but varying lattice constant a e {1 /2, 1 /3} demonstrating the 
uv convergence. The sizes of the dvr basis sets are Lk c /n = 30, 
40, and 50 (blue pluses) and 60, and 90 (red circles) respectively. 
For the blue pluses, the corresponding number of harmonic 
oscillator wave functions suggested in Refs. [i 2] (see also 
Eqs. Q), would be N = Lk c /4 = Ln/4a w 24, 31, 39 and 47, 71 
respectively. Notice that the size of the dvr basis set can be 
reduced by factor of 2 to Lk c /2n = 15, 20, 25 (blue) and 30, 45 
(red) respectively, by imposing Dirichlet boundary conditions, 
however, in that case, states not localized to a single cell will 
not be reproduced. 



In Fig. [i]we show the energy differences between the 
eigenvalues of the Hamiltonian (25 1 and hcu(n+ 1/2). 
These "errors" are indicative only of the energy shifts 
due to the tunneling between neighbouring cells in the 
case of periodic boundary conditions, as one can judge 
by comparing systems with different lengths at the same 
energy, when the tunneling matrix elements are similar. 
The results for the lowest 2/3 of the spectrum are, for 
all practical purposes, converged in the dvr method, 
and the harmonic oscillator basis set is worse in this 
case. With N = Lk c /4 w 24 one can obtain at most 10 
states or so with a reasonable accuracy in this reduced 
interval on the x-axis with periodic or Dirichlet bound- 
ary conditions, if one were to follow the prescription of 

Refs. mm. 

In Fig. [2] we demonstrate the uv and ir exponen- 
tial convergence of the dvr method for an asymmetric 
short-range potential with analytic wavefunctions. Note 
that both ir and uv errors scale exponentially until ma- 
chine precision is achieved - AE oc exp(— 2k OQ L) (left) 
and AE oc exp(— 2k c ro) (right) respectively, where ro is 
potential dependent and koo is determined by the bound 
state energy E = — h 2 k^/2m. These exponential scal- 
ings follow from simple Fourier analysis (uv) and band 
structure theory (ir) for short-ranged smooth potentials. 
Note in particular that the linear uv scaling (right) dif- 



Figure 2. (color online) Here we demonstrate the exponen- 
tial convergence of the periodic dvr basis for the energy of 
the bound states of analytically solvable Scarf II potential 
V(x) = [a + b sinhx]/ cosh 2 x. With a = 7/2 and b = -1 1 /2, the 
potential has three bound states - E n = —(3 — n) 2 /2 (shown 
in black, blue, and green from left to right respectively). In 
the left plot, we demonstrate the ir convergence for increasing 
L with fixed k c , while in the right plot we demonstrate the 
uv convergences for increasing k c for fixed L. The various 
values for k c e {5,10,15,20} (left) and L e {5,15,25,35} (right) 
correspond to dotted, dot-dashed, dashed, and solid lines with 
increasing convergence respectively. 



fers from the quadratic empirical dependence discussed 
in |UJ. 

In summary, the dvr basis seems ideal for nuclear 
structure calculations. It is near optimal in size, and 
can deliver results with exponential convergence. This 
dvr basis shares the important advantages of the ho 
basis set: efficiently separating out the center-of-mass 
motion using Jacobi coordinates (with the added benefit 
of not needing to evaluate Talmi-Moshinsky coefficients), 
utilizing symmetries to reduce the basis size (spherical 
with the Bessel function dvr), and coverage of the low- 
energy ho spectrum. Moreover, matrix elements are 
easy to evaluate - the potential matrix is diagonal for 
local potentials (no overlap integrals are needed - see for 
example Eq. (19)), and the kinetic energy matrix is sparse 
and explicitly expressed analytically. Furthermore, the 
uv and ir convergence properties of the basis appear 
on a equal footing, and are clearly expressed in terms 
of the momentum-space projection and finite box size, 
allowing for simplified and sound convergence analysis. 

We thank G.F Bertsch for discussions. Ideas de- 
scribed here were developed with support under us 
doe grants DE-FG02-97ER41014, DE-FC02-07ER41457, de- 
FG02-00ER41132, and DE-FG02-00ER41132. 
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